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A  NONLINEAR  PROGRAMMING  ALGORITHM  FOR  OPTIMIZING 
CONVENTIONAL  BEAMFORMER  SHADING  WEIGHTS 


1.  INTRODUCTION 


Environmental  adaptation  of  sonar  systems  is  the  process  of  improving  the  performance  of 
a  sonar  system  given  some  known  environmental  conditions  (such  as  measured  environmental 
noise).  Since  a  well-designed  system  is  usually  environmental-noise  limited,  the  use  of  this 
knowledge  of  the  environment  has  the  potential  for  large  performance  improvements.  An  effort 
has  been  undertaken  to  develop  an  optimization  system  to  improve  the  performance  of 
conventional  beamformers  under  known  environmental  noise  conditions.  The  resulting 
optimization  algorithm  shows  great  potential  for  improving  passive  sonar  system  performance 
against  unknown  signals  when  the  (element-level)  noise  performance  is  known.  This  algorithm 
is  directly  applicable  to  the  conventional  delay-and-sum  beamformer  and  thus  requires  minimal 
system  impact  for  implementation. 

The  use  of  delay-and-sum  beamformers  is  commonplace  in  passive  sonar  systems  design. 
These  are  often  referred  to  as  conventional  beamformers.  The  conventional  beamformer 
functions  by  applying  a  time  delay  to  each  sensor  and  multiplying  each  of  these  by  a  prescribed 
weight  before  summing  the  element  responses  to  provide  the  beamformer  output.  Traditionally, 
the  beamformer  is  designed  so  that  the  time  delays  maximize  the  sensitivity  of  the  array  to  a 
desired  (look)  direction,  thus  achieving  robust  signal  gain;  and  the  corresponding  weightings  are 
chosen  to  provide  low  sensitivity  to  energy  that  does  not  match  with  the  incident  wave  (low 
sidelobes),  which  achieves  low  noise  response. 

In  the  absence  of  specific  knowledge  of  the  noise  field  encountered,  the  conventional 
beamformer  design  strategy  provides  robust  performance.  However,  components  of  the  noise 
field  are  often  predictable,  particularly  in  the  case  of  arrays  that  employ  long-duration  time 
averaging.  Those  arrays  have  the  opportunity  to  observe  the  noise  field  (as  seen  by  the  array 
elements)  in  real  time.  In  this  case,  adaptive  beamforming  techniques  are  often  employed, 
wherein  the  beamformer  converts  each  of  the  received  sensor  outputs  into  a  digital  form  and 
employs  frequency-domain  calculations  to  cancel  those  components  of  the  received  signal  that 
appear  to  be  attributable  to  noise.  Such  calculations  provide  a  dramatic  change  of  the  structure 
of  the  beamformer  to  provide  this  noise  removal.  As  an  alternative  means  of  reducing  the  impact 
of  noise  on  the  beamformer,  environmentally  adaptive  methods  of  conventional  beamforming  are 
available,  where  the  conventional  beamformer  structure  is  maintained,  and  the  shading  weights 
found  within  that  structure  are  optimized  for  best  performance  against  the  noise  field. 
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2.  OBJECTIVE  FUNCTION  FOR  CONVENTIONAL 
SHADING  WEIGHT  OPTIMIZATION 


The  problem  of  conventional  beamformer  shading  weight  optimization  has  been  recently 
formulated  from  a  deflection  optimization  perspective  (see  reference  1).  Deflection  is  a 
measurement  of  the  statistical  variation  between  two  hypotheses  in  a  binary  (two-state)  statistical 
hypothesis  test.  In  the  case  of  passive  arrays,  these  states  are  the  state  where  signal  is  present 
and  the  state  where  signal  is  absent.  The  deflection  coefficient  is  a  parameter  that  provides  a 
measurement  of  this  variation.  For  the  case  of  distinguishing  signal-plus-noise  (hypothesis  1) 
from  noise  alone  (hypothesis  0),  one  assumes  that  both  processes  are  governed  by  Gaussian 
random  variables  with  the  same  variance.  Such  an  assumption  is  reasonable  in  the  case  of  small 
signal-to-noise  ratio  (SNR)  that  is  seen  in  passive  sonar  system  applications.  In  that  case,  the 
deflection  coefficient  is  given  by 

d  _  Msn  ~  Mn  t  (1) 

° N 


where  /Jsn  is  the  mean  under  the  signal-plus-noise  hypothesis,  /Un  is  the  mean  under  the  noise- 
only  hypothesis,  and  aN  is  the  standard  deviation  of  the  noise  output  (recall  that  this  is  the  same 
as  that  for  signal-plus-noise).  In  a  standard  square-law  detector,  the  standard  deviation  of  the 
noise-only  hypothesis  crN  is  given  by  (see  reference  2) 

crN  =  ^7  ^HQ(<oj[V%(e>)dca ,  (2) 

where  Hq(co)  is  the  pre-detection  filter  transfer  function,  VN(co)  is  the  beamformer  output  noise 
spectrum,  and  T  is  the  averaging  time.  For  the  sequel,  the  pre-detection  filter  transfer  function  is 
assumed  to  be  unity;  that  is,  H0((o)  =  1 .  Such  an  assumption  is  reasonable  in  the  absence  of 

specific  knowledge  of  the  source  signal.  Furthermore,  the  averaging  time  is  assumed  to  be  large 
compared  to  the  correlation  time  of  the  beamformer  output  noise.  The  mean  value  of  the  noise- 
only  hypothesis  is  given  by  (see  reference  2) 

2 

*£>.(«!  VN(a)da,  (3) 

2n  ^ 

with  a  similar  expression  holding  for  the  signal-plus-noise  mean  value.  The  statistical 
independence  of  signal  and  noise  implies  that  Vs(<y)=  Vsn((0)  -  Vn(co),  so  that  the  deflection 
coefficient  (equation  (1))  reduces  (in  the  case  of  a  unity  pre-detection  filter)  to  the  following: 


(4) 


3 


The  goal  of  the  passive  sonar  detection  problem  is  now  restated  as  a  maximization  of  the 
deflection  coefficient  (equation  (4))  under  the  signal  and  noise  conditions  of  interest. 


In  the  passive  sonar  application,  the  signal  is  unknown,  yet  the  array  is  steered  to  a  prescribed 
direction  6.  Thus,  the  source  is  modeled  as  a  plane  wave  arriving  from  direction  (9  with 
unknown  spectrum  S(a>).  Consider  a  passive  array  of  M  sensors  with  a  conventional  delay-and- 
sum  beamformer.  Let  ks  be  the  corresponding  wavevector  of  the  plane  wave  from  direction  6 . 
The  beamformer  output  power  spectrum  in  the  direction  <9is  represented  by  the  autocorrelation 
of  the  time  output  of  the  delay-and-sum  beamformer.  Specifically,  this  is  given  by 

V(a>)  =  £*E[v(0  v{t  +  t)]  e~ia)T  dr ,  (5) 

where  the  beamformer  output  voltage  v(t)  is  given  by 

M 

v(0=Zwffl  Um(t-Tm),  (6) 

m=\ 


with  um(t)  representing  the  time -domain  output  of  the  m-th  sensor  and  rm  representing  the 
appropriate  time  delay  to  match  the  wavevector  k$. 

For  the  case  of  sensors  receiving  noise-only  input,  the  beamformer  output  power  spectrum  in 
equation  (5)  has  been  shown  (reference  3)  to  reduce  to 

^(^)  =  WruV)M(<u)U(®)W,  (7) 

where  superscript  T  represents  transpose,  (•)*  represents  the  complex  conjugate  transpose,  and  W 
and  U(a>)  are  given  by  the  following: 

W  =  [w1,w2,---,wA/]r,  (8) 


and 


U(<y)  =  Diag{exp(/ks  •x]),exp(/k5  •x2),---,exp(/k5  -xM)}.  (9) 

The  matrix  M(<w)  is  a  matrix  of  sensor  noise  measurement  cross-correlations;  that  is,  the  0V)-th 
component  of  M(<y)  is  given  by  the  cross-correlation  of  the  noise  output  of  sensor  i  with  sensor  j. 
For  the  case  of  sensors  receiving  signal-only  input  with  spectrum  S(a>)  in  the  form  of  a  plane 
wave  in  direction  6,  the  beamformer  output  power  spectrum  in  equation  (5)  has  been  shown 
(reference  4)  to  reduce  to 
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(  M 


Vs(co)  =  S(co)G 


^2 


2>„(*) 


wn 


\m=\ 


(10) 


where  G  is  the  power  gain  of  an  individual  sensor  and  dm{9 )  is  the  directivity  loss  of  the  m-th 
sensor  when  receiving  a  signal  in  the  direction  9.  For  example,  omnidirectional  sensors  have 
dm(9 )  =  1 ,  and  cosine-directive  velocity  sensors  have 

ffks((9)-n,„)/*o,  for  ks(6)  nm  >0 

t  .  •  (U 
|  0,  otherwise 

where  n,„  is  the  sensor  outward  normal  and  k0  is  the  free-space  wavenumber.  Equations  (7)  and 
(10)  provide  the  power  spectra  required  for  general  evaluation  of  the  deflection  coefficient 
expressed  in  equation  (4).  Together,  these  equations  represent  the  analytical  form  of  the 
conventional  passive  sonar  array  performance. 

The  goal  of  conventional  beamformer  shading  weight  optimization  is  to  maximize  the 
performance  (as  measured  by  the  deflection)  through  the  adjustment  of  the  shading  weights  wm. 
This  is  accomplished  by  minimizing  the  denominator  of  the  deflection  while  holding  the  value  of 
the  numerator  fixed.  This  problem  is  stated  in  the  form  of  a  nonlinear  optimization  problem  as 
follows: 


min  f  Vl(co)dco 

W  J-oo 


such  that  (12) 

M 

m= 1 

The  constraint  expression  in  (12)  is  a  direct  result  of  holding  the  beamformed  signal  power 
output  expression  (10)  constant  with  respect  to  the  available  parameters  wm.  Note  that  the  sensor 
gain  G  and  the  signal  spectrum  S(a>)  are  not  available  as  degrees  of  freedom  in  the  optimization 
and,  therefore,  are  taken  as  unknown  constant  parameters  with  respect  to  the  optimization 
problem  (12).  The  optimization  problem  (12)  is  similar  in  intent  to  the  goal  of  Minimum 
Variance  Distortionless  Response  (MVDR)  beamforming,  but  the  optimization  problem  varies  in 
two  important  respects.  First,  the  weights  are  constrained  to  be  real  and,  second,  a  single  set  of 
weights  is  found  over  the  entire  frequency  band  of  interest.  These  differences  make  the 
nonlinear  programming  approach  amenable  to  the  conventional  beamformer  design  and  thus  it  is 
not  an  “adaptive  beamformer”  approach,  although  the  resulting  weights  are  adaptively 
determined.  In  practice,  the  limits  of  integration  in  (12)  are  restricted  to  the  expected  frequency 
band  of  interest,  thus  providing  a  limiting  case  (as  bandwidth  goes  to  zero)  of  the  frequency 
domain  beamformer  optimization  problem  that  is  solved  in  MVDR. 
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3.  SOLVING  THE  OPTIMIZATION  PROBLEM  USING 
EXISTING  NUMERICAL  TECHNIQUES 


The  numerical  optimization  of  nonlinear  objective  functions  under  constraints  (the 
nonlinear  programming  problem)  is  a  well-studied  problem  with  many  available  solution 
techniques.  Unfortunately,  such  problems  are  often  prone  to  multiple  local  minima,  and 
numerical  methods  that  do  not  start  near  the  solution  may  converge  to  one  of  these  locally 
optimal  solutions  with  no  warning.  For  the  special  case  of  nonlinear  programming  problems  that 
are  convex  optimization  problems,  all  local  minima  are  guaranteed  to  be  global  minima ,  so  the 
problem  of  local  solutions  no  longer  exists.  For  these  problems,  the  quadratic  programming 
(QP)  and  sequential  quadratic  programming  (SQP)  methods  are  often  utilized  due  to  their  rapid 
convergence  properties  and  ease  of  implementation.  QP  methods  are  applicable  only  to  a  further 
subset  of  problems  in  which  the  objective  is  quadratic,  whereas  SQP  methods  apply  a  process  of 
sequentially  approximating  the  problem  as  a  quadratic  objective  in  a  convergent  manner.  In  this 
section,  the  utility  of  conventional  SQP  methods  in  solving  the  optimization  problem  (equation 
(12))  is  examined.  A  thorough  overview  of  SQP  methods  is  given  in  reference  5. 

The  objective  function  and  simple  linear  constraints  shown  in  equation  (12)  are  a  convex 
optimization  problem.  By  definition,  an  optimization  problem  is  convex  if  the  objective  function 
is  a  convex  function  and  the  set  of  feasible  points  (the  set  of  variables  that  meet  the  constraints) 
is  convex.  For  the  problem  under  consideration  (i.e.,  equation  (12)),  the  objective  is  convex 
since  it  is  the  square  of  a  quadratic  form  (equation  (7)  shows  that  Vjy ((d)  is  a  quadratic  form  in  the 
weights  W),  and  the  constraint  set  is  a  convex  linear  combination  of  the  variables.  Thus,  problem 
(12)  is  a  convex  nonlinear  programming  problem  and  is  thus  suitable  for  solution  using  SQP 
methods. 

Consider  a  general  nonlinear  programming  problem  as  an  objective  function  and 
constraints  as  follows: 

min/(x) 


such  that  (13) 

C(jc)  <  0 . 

Any  general  nonlinear  programming  problem  can  be  placed  in  this  form  by  incorporating  all 
equality  constraints  into  the  objective  function.  The  Lagrangian  function  of  the  optimization 
problem  (13)  combines  the  objective  and  constraints  into  a  single  function  L(x,X),  where  X 
represents  the  Lagrange  multipliers  associated  with  (13).  The  general  SQP  method  considers  a 
sequence  of  estimates  of  x  (given  by  xk),  where  at  each  estimate  a  Taylor  series  expansion  of  the 
objective  is  used  to  obtain  the  quadratic  function 

f(xM)  *  f(xk)  +  Vf{xk)Tdk+\dTk  V2L(xk,Ak)dk,  (14) 
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where 


dk 


~xk+ 1  ~xk 


(15) 


represents  the  step  taken  to  get  to  the  next  iterate.  It  can  be  shown  (see  reference  6)  that  this 
iteration  scheme,  when  convergent,  leads  to  a  solution  that  satisfies  the  first  and  second  order 
conditions  for  optimality  (the  Karush-Kuhn-Tucker  conditions).  The  convergence  is  guaranteed 
whenever  the  Hessian  of  the  Lagrangian  maintains  positive-definiteness.  The  numerical 
minimization  of  quadratic  equation  (14)  is  easily  performed  using  standard  QP  solution 
techniques  (see  reference  7). 

Most  commercial  general-purpose  optimization  codes  for  nonlinear  programming  solve 
equation  (14)  iteratively.  The  computational  savings  is  that  a  single  algorithm  easily  handles 
general  objective  functions  and  constraints.  This  makes  it  appropriate  for  general-purpose  codes. 
However,  the  user  must  still  define  the  objective  function,  the  gradient  function,  and  the  Hessian 
in  a  manner  that  is  suitable  for  the  general-purpose  solver.  One  of  the  most  reliable  of  the 
general-purpose  SQP  solvers  is  NPSOL  (reference  8)  from  Stanford  University’s  Systems 
Optimization  Laboratory.  NPSOL  solves  equation  (14)  for  general  objective  functions  and 
general  constraints.  The  user  must  provide  a  separate  subroutine  to  evaluate  the  objective 
function  and  the  gradient  at  an  arbitrary  iteration  step.  The  Hessian  of  the  Lagrangian  is 
evaluated  approximately  using  a  Broyden  Fletcher  Goldfarb  Shanno  (BFGS)  update  formula  (see 
reference  9).  Such  an  update  provides  a  computationally  reliable  technique  to  improve  estimates 
of  the  Hessian. 

A  numerical  implementation  of  array  shading  weight  optimization  using  NPSOL  was 
employed  as  the  first  component  of  this  effort.  The  results  of  this  algorithm  are  found  in 
references  4  and  10.  For  the  problem  under  consideration  (equation  (12)),  the  required  objective 
function  and  gradient  function  are  written  as 

/(W)  =  0\VrU*(cy)M(a>)U(«)w)2^,  (16) 

and 

V/(W)  =  4  r>l(wrU*(<y)M(<y)U(<sj)w)u*(fy)M(<y)U(iy)W dco ,  (17) 


where  the  integrals  are  evaluated  using  the  method  of  overlapping  parabolas,  and  the  frequency 
limits  of  [a*,  co\]  represent  the  limiting  band  of  the  sonar  system  design.  The  numerical 
implementation  of  this  objective  within  NPSOL  worked  well  for  small  arrays.  As  the  number  of 
elements  (size  of  W)  grew,  the  performance  of  NPSOL  degraded.  In  particular,  the  ability  of 
NPSOL’ s  BFGS  algorithm  to  handle  the  updating  of  such  a  large  Hessian  (the  Hessian  size  is 
M  x  M  for  an  M-element  array  problem)  faulted.  In  particular,  for  arrays  larger  than  100 
elements,  the  Hessian  update  algorithm  actually  diverged,  producing  a  total  lack  of  convergence. 
Numerical  investigations  showed  that  the  problem  is  inherent  in  using  a  BFGS  update  of  the 
Hessian,  and  a  method  without  such  a  requirement  was  sought. 
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An  alternative  approach  to  the  use  of  a  general  solver  for  the  solution  of  the  QP 
subproblem  (14)  is  to  develop  a  specialized  algorithm  that  lacks  the  flexibility  to  handle  general 
objectives  but  is  tuned  so  that  the  Hessian  is  included  exactly.  The  main  difficulty  with 
including  the  Hessian  directly  is  that  the  Lagrangian  is  usually  much  more  complex  than  the 
simple  objective.  By  limiting  attention  to  the  array  shading  problem  (12),  the  constraints  are  all 
linear;  hence,  the  Hessian  of  the  objective  is  identical  to  the  Hessian  of  the  Lagrangian  (all 
second  derivatives  of  constraint  terms  are  necessarily  zero).  Thus,  for  the  specific  array  shading 
optimization  problem,  the  Hessian  is  explicitly  available  as 


V2/(W)  =  4£'U*(<y)M(iy)U(6>)(wrU*(iy)M(<y)U(6;)w)j<y 

+  8  r'U’(6))M(6>)U(6>)WWrU(to)M(6))U*((y)^, 


(18) 


which  is  readily  computed  at  each  iteration.  A  general-purpose  nonlinear  programming  solver 
does  not  allow  such  details  to  be  input  by  the  user.  By  developing  a  specialized  code  that  solves 
only  the  array  shading  weight  optimization  problem,  the  numerical  issues  associated  with 
approximation  updates  to  the  Hessian  will  no  longer  be  a  concern,  since  the  Hessian  will  be 
evaluated  directly  at  each  step.  It  is  noted  that  such  an  approach  is  not  advisable  for  small 
problems,  since  the  evaluation  of  (18)  has  considerable  computational  costs  relative  to  applying 
the  computationally  efficient  BFGS  method. 


4.  ARRAY  SHADING-SEQUENTIAL  QUADRATIC  PROGRAMMING 

(AS-SQP)  ALGORITHM 


4.1  MATHEMATICAL  DESCRIPTION 

The  Array  Shading-Sequential  Quadratic  Programming  (AS-SQP)  algorithm  was 
developed  to  provide  robust  numerical  solutions  to  the  specific  array  shading  optimization 
problem  (12)  in  cases  where  general-purpose  solvers  fail  to  converge.  The  mathematical 
derivation  of  AS-SQP  is  developed  for  a  general  convex  objective  function  with  a  single  convex 
linear  equality  constraint  and  simple  bound  constraints  on  the  variables  (array  weights). 

Application  to  the  problem  (12)  is  accomplished  by  evaluating  expressions  (16),  (17),  and  (18)  to 
the  appropriate  components  of  the  algorithm.  Consider  the  following  nonlinear  optimization 
problem: 

min  F( w)  (19a) 

w  gRm 


subject  to  the  linear  equality  constraint 

M 

2X=i  (19b> 

m= 1 
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and  the  lower-bound  and  upper-bound  constraints  on  the  array  weights 


0  <  wm  <  1,  m  =  1,  2, ...,  M.  (19c) 

In  the  above  formulas,  the  array  weights  are  expressed  in  terms  of  the  vector  w  =  [w„  w2,  ..., 
wm]t.  This  formulation  is  consistent  with  the  nonlinear  programming  problem  formulated  in  the 
preceding  section.  Note  that  non-unity  element  contributions  to  the  steer  direction  (i.e., 
directivity  effects)  can  be  accounted  for  by  adjusting  the  definition  of  “weight”  to  include  the 
conventional  weight  combined  with  the  directive  component  in  the  steered  direction  (i.e., 
wm  <-  dm(0)wm ).  For  simplicity  of  exposition,  the  case  of  unity  directivity  (omnidirectional 
sensors)  is  considered  in  this  mathematical  derivation;  however,  the  computational 
implementation  includes  the  directivity  effects. 

Generally,  an  SQP  algorithm  involves  two  parts:  the  determination  of  the  search  direction 
8  obtained  from  the  solution  of  a  QP  subproblem  and  the  determination  of  the  step  size  along  the 
search  direction  8.  At  the  end  of  each  iteration,  the  weight  vector  w  is  updated  and  the  procedure 
is  repeated  until  either  a  convergent  result  is  obtained  or  the  maximum  number  of  allowed 
iterations  is  reached.  At  each  iteration  of  the  SQP  algorithm,  the  gradient  g  and  the  Hessian  H  of 
the  objective  function  F( w)  are  evaluated  at  the  current  iterate  w(k\  where  k  is  the  iteration  index. 
The  search  direction  8(A)  at  the  point  yv{k)  is  then  determined  through  the  solution  of  a  QP 
subproblem.  The  subproblem  involves  the  minimization  of  the  quadratic  function  approximation 
gW  to  the  objective  function  at  the  point  w(k)  subject  to  constraints  (19b)  and  (19c). 

The  AS-SQP  algorithm  is  explicitly  focused  on  the  problem  of  optimal  shading  weights. 
Thus,  rather  than  a  generic  set  of  constraints,  only  linear  shading  weights  of  the  form  (19b)  are 
considered  (along  with  bound  constraints  on  the  individual  weights).  The  QP  subproblem  for 
AS-SQP  is  given  as 

min  Qik)(b)  (20a) 

beRM 

subject  to  the  linear  equality  constraint 

IX  =o  <20b> 

m= 1 


and  the  lower-  and  upper-bound  constraints  on  the  components  of  the  search  direction  vector 

-1<  Sm<\,  m=l,2,...,M.  (20c) 

In  expression  (20a),  the  objective  function  is  defined  as 


Q{k\b)  =  ^8rH(t)8  +  dTgik)  +  F( w{k) ), 


(21) 


where  H(i)  and  g(A)  denote  the  Hessian  and  gradient  of  F( w),  respectively,  evaluated  at  the 
current  iterate  w(k\  Constraints  (20b)  and  (20c)  for  the  search  direction  8  are  consistent  with 
array  weight  constraints  (19b)  and  (19c). 

To  derive  constraints  (20b)  and  (20c)  in  the  QP  subproblem,  the  search  direction  vector  is 
expressed  as 

5  =  w-w^.  (22) 

The  components  of  8  are  given  as 

(23) 

If  the  M  components  of  equation  (23)  are  summed,  linear  constraint  (20b)  is  obtained;  i.e., 

M  M  M 

<24> 

m-\  m= 1  m-1 

where  expression  (19b)  has  been  applied  to  the  sums  involving  w  and  Lower-  and  upper- 
bound  constraints  (20c)  of  the  components  of  the  search  direction  vector  can  be  derived  from 
formulas  (19c)  and  (23).  From  constraint  (19c),  one  has 

-1<  wm  <1,  m  =  1,2,. ..,M.  (25) 

Substitution  of  expression  (23)  into  equation  (25)  yields  the  variable  constraints  (20c)  in  the  QP 
subproblem. 

Once  the  solution  to  the  QP  subproblem  (20)  is  obtained,  the  new  iterate  w(W)  is  given  by 

w(*+1  )=w(*)+a(*)g(*)  (26) 

where  the  non-negative  scalar  ot®  is  the  step  size  to  be  determined.  In  the  QP  subproblem  (20), 
because  the  quadratic  function  approximation  to  F(w)  at  w  =  is  minimized  at  w  =  +  8^, 

this  minimum  corresponds  to  a  step  size  of  ol®  =  1  in  expression  (26).  However,  because  the 
components  of  w^+1)  must  lie  within  the  constraint  region  (19c),  limitations  are  placed  on  dk\ 

The  second  part  of  the  SQP  algorithm  involves  the  determination  of  the  step  size  of®  in 
expression  (26).  The  components  of  the  new  iterate  w(*+|)  must  satisfy  constraints  (19c);  i.e., 

0  <  w(k)  +  a(k)S(mk)  <  1,  i»  =  l,2,...,M.  (27) 

In  addition,  the  step  size  must  satisfy 
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0  <a(k)  <1. 


(28) 


Inequality  (27)  places  constraints  on  the  step  size  ct®  in  order  to  ensure  that  the  components  of 
w(*+1)  stay  within  the  constraint  space.  Because  the  new  iterate  w^+1^  must  satisfy  the  equality 
constraint  (19b),  the  upper  bound  of  inequality  (27)  will  never  be  reached.  The  mth  component 
of  the  argument  in  expression  (27)  will  reach  or  exceed  its  lower-bound  constraint  for  the 
following  step  sizes: 


„(*) 


a{k)>-^frr,  S(mk)>0,  m  =  \,2,...,M , 

8{k) 


„(*) 


a{k)<-^rr,  S(mk)<0,  m  =  . 

sly 


(29) 


(30) 


Condition  (29)  is  trivial  because  cfk)  is  already  restricted  to  the  interval  [0,1].  Therefore,  the  step 
size  dk)  in  the  optimization  problem  (27)  is  determined  from  the  components  of  w(*+1)  with 
negative  search  direction  components.  Define  the  set  I \  as  follows: 


O)-  <31> 

Therefore,  from  formulas  (30)  and  (31),  the  maximum  allowable  step  size  at  the  kth  iteration  is 
given  by 


a 


(k) 


f 

=  min 

meL 


,(*)' 


(32) 


The  AS-SQP  algorithm  also  contains  a  user-defined  step  size  au,  which  serves  as  the 
default  step  size  during  each  major  iteration,  where  0  <  ccu  «  1.  The  user-defined  step  size  is 
small  because  the  gradient  and  Hessian  of  F(w)  will  vary  as  w  moves  along  the  search  direction. 
At  each  major  iteration,  the  step  size  is  determined  as 


=  min<! 


a,, ,  min 

«e/; 


w, 


(33) 


At  the  end  of  each  iteration,  the  components  of  the  updated  weight  vector  wa+1)  are  checked  to 
see  if  they  are  active  (i.e.,  to  see  if  any  lie  on  the  lower  boundary  of  the  constraint  region  (19c)). 
For  each  active  component  of  w(A+1),  the  lower  bound  of  the  corresponding  search  direction 
component  is  set  to  zero  (i.e.,  0  <  8^^ 1  <1  Vma  ' 1  =  0 ).  The  QP  subproblem  is  then 

recalculated  with  the  updated  weight  vector  and  direction  constraints.  The  above  procedure  is 
repeated  until  either  a  convergent  result  is  obtained  or  the  maximum  number  of  iterations  is 
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reached.  A  complete  flowchart  of  the  AS-SQP  numerical  algorithm  is  presented  in  the  appendix. 
Note  that  the  addition  of  non-unity  directivity  complicates  the  expressions  somewhat;  however, 
the  formulation  is  identical  to  that  described  above. 


4.2  COMPUTATIONAL  IMPLEMENTATION 

The  numerical  implementation  of  the  AS-SQP  algorithm  has  been  performed  as  a  software 
program  called  NumShade  (abbrevation  for  Numerical  Array  Shading  Weight  Optimization). 
NumShade  was  written  in  FORTRAN  95  with  a  Windows-based  graphical  user  interface  (GUI), 
using  Lahey  FORTRAN  95  version  5.7  along  with  Lahey’s  Winteracter  development 
environment.  The  entire  package  (software  with  GUI)  is  available  as  a  single  Windows 
executable  file  that  will  run  on  any  computer  under  the  Windows  operating  system,  provided 
there  is  enough  memory  available  on  the  computer  for  the  calculations.  Dynamic  memory 
allocation  within  FORTRAN  95  was  employed  to  improve  this  portability.  Thus,  computer 
systems  with  little  memory  will  suffice  for  solving  small  array  problems  (few  elements  and/or 
low  bandwidth)  and  a  larger  memory  computer  can  be  employed  to  solve  a  larger  array  shading 
problem.  This  scalability  of  problem  size  is  inherent  in  the  software  structure  and,  thus,  is 
transparent  to  the  user. 

As  part  of  the  computational  implementation  of  the  AS-SQP  algorithm,  expressions  for  the 
objective  function  and  its  gradient  and  Hessian  were  derived  and  evaluated  directly  in  the 
software.  This  eliminates  the  need  for  numerical  methods  to  approximate  the  derivatives  and 
thus  removes  one  of  the  major  drawbacks  of  using  standardized  numerical  software.  These 
evaluations  are  accomplished  in  a  secondary  subroutine  that  evaluates  expressions  (16),  (17),  and 
(18)  at  a  given  iteration  using  the  method  of  overlapping  parabolas  to  provide  the  numerical 
evaluation  of  the  integrals.  The  GUI  of  the  software  is  shown  in  figure  1.  A  second  software 
program  that  employs  NPSOL  (for  use  with  small  array  problem)  is  also  available  with  an 
identical  GUI.  The  INPUTS  side  of  the  GUI  takes  in  filenames  for  element  position  data, 
element  normal  data,  and  noise  cross-correlation  data.  These  filenames  may  be  typed  in,  or  the 
appropriate  search  button  can  be  used  to  open  a  standard  Windows  file-chooser  box.  All  of  the 
input  files  are  ASCII  text  files  that  contain  the  data  specific  to  an  array.  The  Sensor  Locations 
file  consists  of  M  +  1  lines,  the  first  line  containing  the  number  of  elements,  and  the  remaining 
lines  containing  the  (x,  y,  z)  positions  of  the  elements  (separated  by  spaces).  The  Sensor 
Normals  file  is  identically  formatted,  except  that  the  positions  are  replaced  with  the  (x,  y,  z) 
components  of  the  unit  outward  normal  to  each  sensor.  This  information  is  used  to  determine  the 
sensor  directivity  response.  The  default  version  of  the  software  assumes  cosine-directive 
elements;  other  element  directivity  patterns  are  readily  incorporated  in  special  releases  of  the 
software.  The  Noise  Cross-Corr  file  contains  the  element-to-element  noise  cross-correlation  for 
each  frequency  of  interest.  This  file  is  by  far  the  largest  file  needed  to  run  the  program.  The 
formatting  of  the  file  is  as  follows:  the  first  line  is  frequency  1  value  and  the  next  line  is  2 *M 
numbers  representing  the  frequency  1  component  of  the  complex  noise  cross-correlation  of 
sensor  1  with  each  other  sensor,  in  order.  The  complex  components  are  represented  in  pairs  of 
real  component  followed  by  imaginary  component,  with  a  space  between  the  two,  then  another 
space  before  the  next  complex  pair.  The  successive  lines  continue  this  cross-correlation  for  each 
of  the  sensors  until  frequency  1  is  exhausted;  then,  the  process  is  repeated  for  each  frequency. 
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The  resulting  file  contains  F*(M+ 1)  rows,  where  F  is  the  number  of  frequencies  and  M  is  the 
number  of  elements. 


JU 


-  INPUTS  - 


SEARCH 


Sensor  Locations  lgeom.dat 
Filename 


Sensor  Normals  Lorm.dat 
Filename 


Ji 


J 


N  oise  Cross-Corr  IxcorLfull^daT 
Filename  1 


First  Element:  |  18l]|  Last  Element:  |  26Q| 

N umber  of  Frequencies:  f~ To| 
r""1^ 7: 

Steering  Angle: 


9011  degs  azimuthal 


elevation 


OUTPUTS 


Shading  Weights 
Filename 

Numerical  Summary 
Filename 


|weight$_optout 

_JJ 

|num_results.out 

3 

.  — ^ACflONS  — 1 

rff|  HELP  I 

PLOT  |  . EXif  | 


Figure  1.  Graphical  User  Interface  for  the  Optimization  Software 


The  user  has  the  option  of  selecting  a  subset  of  the  sensor  elements  (or  subarray)  by 
specifying  a  first  and  last  element.  This  allows  the  entire  database  of  noise  information  to  be 
stored  in  a  single  ASCII  text  file,  and  the  software  accesses  the  components  necessary  for  the 
optimization  that  is  desired.  Furthermore,  the  number  of  frequencies  can  be  specified  as  less 
than  the  number  available.  At  this  time,  the  frequency  list  always  begins  with  the  first  frequency 
in  the  noise  cross-correlation  file.  The  final  setup  input  is  the  array  steering  angle,  which 
assumes  a  traditional  submarine  orientation  with  (azimuthal,  elevation)  pairs  of  (0°,0°)  pointing 
along  the  negative  x-axis,  (90°, 0°)  pointing  along  the  y-axis,  and  (0°,90°)  pointing  along  the  z- 
axis.  A  final  user  input  is  the  selection  of  output  filenames.  The  Shading  Weights  file  contains 
two  columns — element  number  and  final  shading  weight.  The  Numerical  Summary  file  contains 
summaries  of  the  numerical  performance  of  the  algorithm,  which  is  used  by  advanced  users  for 
assessing  convergence  performance. 

Once  the  files  are  set  up  and  the  appropriate  fields  are  all  filled  in  on  the  GUI,  the  user 
presses  the  RUN  button.  After  the  data  files  are  loaded  into  memory,  a  secondary  warning 
window  pops  up  to  inform  the  user  that  the  run  is  set  up.  When  the  OK  button  on  that  window  is 
pressed,  the  optimization  calculation  begins.  The  output  data  files  can  be  read  in  to  any  offline 
program  of  the  user’s  choosing.  There  is  a  very  basic  plotting  facility  under  the  PLOT  button 
that  will  display  a  window  with  the  nominal  (unity  shading)  objective  and  the  resulting 
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optimal  objective  plotted  over  the  frequency  band  of  the  calculation.  It  is  not  recommended  that 
this  plotting  be  used  for  anything  other  than  “quick-looks”  at  the  results  to  determine  that  the 
program  is  functioning  properly.  The  HELP  button  contains  contact  information  for  the  author. 
The  EXIT  button  closes  the  program  window. 


5.  NUMERICAL  RESULTS 


A  set  of  test  data  from  the  WHITEFISH  structural  acoustic  model  was  taken  under  a 
conformal  velocity  sonar  array  test  in  1 996  at  the  Intermediate  Scale  Measurement  System 
(ISMS).  This  data  set  was  collected  to  address  many  issues,  including  a  high  spatial  resolution 
measurement  of  the  structural  noise  field  on  a  large-scale  array  model.  For  the  purposes  of  the 
array  shading  optimization  algorithm,  the  data  set  provides  a  convenient  mechanism  for  testing 
the  array  shading  algorithms  on  a  large  array  with  a  reasonable  level  of  complexity  in  the  noise 
field.  The  test  is  described  in  detail  in  reference  10  and  the  reader  is  referred  there  for  specific 
details  of  the  experimental  setup.  The  test  setup  of  the  WHITEFISH  model  for  array  data 
collection  is  shown  in  figure  2.  The  model  contains  a  series  of  both  small  and  deep  frames 
within  a  steel  pressure  hull.  A  sequence  of  25  shakers  were  employed  along  the  model  to 
provide  very  detailed  structural  acoustic  forcing  patterns  over  a  nondimensional  frequency  range 
of  0.5  <  kQa  <  30 ,  where  ko  is  the  wavenumber  in  the  fluid  and  a  is  the  nominal  radius. 


Figure  2.  Array  Geometry  for  the  Test  Data  Set 


The  WHITEFISH  model  hull  was  coated  with  a  pressure-release  material,  and  conformal 
velocity  sensors  were  placed  in  two  array  configurations  as  shown  in  figure  2.  The  data  used  to 
test  the  optimization  algorithms  come  from  the  primary  array  region,  which  consists  of  a  patch 
of  480  regularly-spaced  elements  as  shown  in  figure  3.  The  elements  that  are  colored  red  were 
considered  faulty  elements  and  were  removed  from  the  analysis.  The  primary  area  of  interest  in 
this  numerical  study  is  subarrays  of  varying  sizes  beginning  at  element  181  and  continuing  until 
element  380.  This  provided  arrays  of  sizes  1  x  20  up  to  10  x  20  for  analysis  purposes.  Further- 
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more,  to  reduce  the  size  of  data  files,  only  a  subset  of  the  frequency  band  with  bandwidth  of 
approximately  k0a  =  3.2  was  used  in  the  numerical  computations.  Examination  of  the  data 

showed  that  this  limited  frequency  band  was  sufficient  to  capture  the  complexity  of  the  noise 
field  and  thus  illustrate  the  array  shading  weight  optimization  algorithm. 
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Figure  3.  Map  of  the  Sensor  Locations  Used  in  the  Test  Array 
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Because  NPSOL  was  shown  to  perform  well  on  small  arrays  (see  section  3),  a  subarray  of 
80  elements  (from  numbers  181-260)  was  used  to  compare  the  convergence  of  the  AS-SQP 
algorithm  with  NPSOL.  This  4  x  20  array  provides  enough  aperture  that  noticeable  benefit  is 
expected  from  optimizing  shading  weights  based  on  noise  measurements.  The  performance  of 
each  algorithm  is  considered  by  measuring  the  relative  improvement  in  broadband  SNR  relative 
to  that  achieved  by  a  unity-shaded  array  over  the  same  frequency  band.  The  use  of  SNR  as  a 
performance  metric  is  natural  since  that  is  a  design  goal  of  any  array.  For  the  passive  array,  the 
signal  is  unknown;  however,  since  the  optimization  algorithm  constrains  the  signal  gain  to  be 
constant  (see  equation  (12)),  the  resulting  ratio  of  SNRs  (or  SNR  improvement)  is  independent 
of  the  signal  level.  Thus,  the  improvement  in  SNR  relative  to  unity  shading  is  the  primary 
performance  metric  for  the  shading  algorithms.  The  result  of  the  convergence  versus  iteration 
number  for  the  two  algorithms  is  shown  in  figure  4.  For  this  computation,  a  frequency 
bandwidth  of  k^a  =  2.5  was  used. 


Interation  Number 


Figure  4.  Algorithm  Convergence  Measured  by  SNR  Improvement  Over  Unity  Shading 


From  figure  4,  it  is  clear  that  both  algorithms  converge  to  a  similar  solution  (and  converge  to 
the  same  amount  of  improvement  over  unity  shading).  Both  algorithms  converge  uniformly 
(they  always  improve)  as  is  expected  for  any  SQP  algorithm  applied  to  a  convex  nonlinear 
programming  problem.  Also,  the  AS-SQP  converges  more  rapidly  initially,  yet  the  NPSOL 
algorithm  converges  at  a  steadier  rate.  This  is  attributed  to  the  Hessian  evaluation  strategies 
employed.  Since  NPSOL  continually  re-estimates  the  Hessian  in  an  iterative  fashion,  it  is 
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limited  to  incremental  performance  improvement  at  each  iteration.  AS-SQP,  on  the  other  hand, 
re-computes  the  Hessian  from  “scratch”  at  each  iteration,  thus  greatly  improving  convergence 
when  there  is  a  direction  of  rapid  change  to  be  followed  to  the  next  iteration.  NPSOL  smoothes 
out  these  rapid  changes  by  only  allowing  gradual  shifts  to  the  Hessian,  thus  providing  slower, 
more  stable  convergence.  While  the  number  of  iterations  implies  that  AS-SQP  converges  more 
rapidly,  each  iteration  is  much  more  expensive,  so  that  AS-SQP  actually  converges  (in  wall¬ 
time)  much  more  slowly  than  NPSOL.  This  is  to  be  expected,  since  the  evaluation  of  the 
Hessian  via  equation  (1 8)  is  the  most  computationally  intensive  part  of  either  algorithm.  Since 
both  algorithms  provide  similar  results  and  run-time  is  not  considered  in  this  study,  AS-SQP  was 
used  in  the  remaining  calculations  due  to  its  ability  to  handle  arbitrarily  large  array  shading 
problems. 

The  tremendous  performance  improvement  shown  in  figure  4  leads  to  the  question  of 
whether  the  bandwidth  was  large  enough  to  consider  this  a  true  broadband  performance  result. 
For  a  broad  enough  frequency  band,  any  complexity  in  the  noise  will  be  seen  by  the  array  as 
white  noise  (due  to  the  frequency  integration)  and,  thus,  the  unity  shading  should  provide  near¬ 
ideal  performance.  To  answer  the  bandwidth  question,  an  examination  of  the  performance 
improvement  relative  to  frequency  bandwidth  was  conducted  using  the  same  4  x  20  subarray. 
The  results  of  both  decreasing  and  increasing  the  bandwidth  are  plotted  in  figure  5.  From  the 
plot,  it  is  clear  that,  although  the  performance  improvement  is  degrading  as  bandwidth  increases, 
the  rate  of  decrease  is  slow  enough  that  substantial  gains  can  be  expected  for  very  large 
bandwidth  problems.  Also,  note  that  the  limit  as  bandwidth  goes  to  zero  approaches  a  value  of 
28  dB,  which  is  the  limit  of  performance  of  the  delay-and-sum  beamformer  for  that  size  of  array 
(superdirectivity  or  non-real  weights  are  the  only  ways  to  improve  performance  beyond  that 
limit).  Thus,  the  array  shading  weight  optimization  software  can  be  employed  to  set 
performance  limits. 
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Figure  5.  Performance  Improvement  as  a  Function  of  Frequency  Bandwidth 


In  addition  to  the  effect  of  performance  degradation  due  to  increased  frequency,  there  is 
also  an  expected  change  in  performance  due  to  changes  in  spatial  extent  (aperture  size).  Since  an 
array  functions  as  a  spatial  filter,  eventually  the  array  size  becomes  large  enough  that  the  noise 
response  appears  spatially  white.  This  phenomenon  is  clearly  illustrated  in  figure  6,  where  a 
variety  of  aperture  sizes  (all  beginning  with  element  number  181)  were  considered.  All  of  the 
calculations  for  this  plot  used  the  koa  =  2.5  bandwidth.  Note  that  the  AS-SQP  algorithm  is 
effective  even  on  arrays  of  200  elements  (twice  the  size  of  the  array  optimization  problem  that 
the  NPSOL  algorithm  handles).  For  large  arrays  (more  than  100  elements),  the  performance 
improvement  degrades  with  array  size,  as  the  spatial  filtering  effectiveness  of  the  aperture 
becomes  a  dominant  mechanism.  This  illustrates  another  area  of  array  analysis  where  the 
optimal  shading  algorithm  is  useful — the  setting  of  array  sizes  to  take  advantage  of  the  spatial 
filtering  effect.  For  the  smaller  apertures  (20, 40,  and  60  elements),  the  array  width  is  so  small 
(all  of  these  have  20-element  columns,  so  the  widths  are  1,  2,  and  3  elements,  respectively)  that 
the  array  shading  cannot  be  used  effectively.  Thus,  there  is  a  minimum  effective  aperture  size 
before  array  shading  optimization  is  merited.  This  size  is  dependent  on  both  the  frequency 
bandwidth  and  the  specific  structure  of  the  noise  profile. 


Number  of  Elements  in  Array 


Figure  6.  Performance  Improvement  as  a  Function  of  Array  Size 


Since  the  array  shading  optimization  algorithms  are  designed  to  minimize  the  noise  response 
Fiv(ft>)  while  simultaneously  maintaining  signal  performance  Vf®),  it  is  expected  that  the 
specific  structure  of  the  noise  field  has  a  large  impact  on  the  resulting  optimal  shading  weights. 
To  examine  the  impact  of  the  specific  structure  of  the  noise  field  on  the  algorithm’s  performance, 
a  comparison  of  optimal  shading  weights  obtained  from  a  simulated  array  with  white  noise 
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(figure  7)  was  computed  and  compared  to  that  from  a  section  of  the  test  array  with  an  identical 
geometry  (figure  8).  In  both  of  these  cases  the  resulting  shading  weights  were  scaled  for 
visualization  so  that  the  maximum  weight  is  at  unity.  The  white  noise  case  in  figure  7  shows  that 
the  array  shading  algorithm  functions  as  expected,  providing  a  standard  tapered  pattern  in  the 
shading  distribution.  The  measured  noise  case  takes  advantage  of  the  element-level  noise 
readings  and  places  higher  weights  on  less-noisy  elements  while  attempting  to  still  distribute  the 
weights  across  the  entire  aperture  to  maintain  array  gain.  This  notion  of  array  gain  is  not  an 
added  feature — the  algorithm  is  only  solving  the  standard  array  shading  optimization  problem  of 
equation  (12). 

To  see  how  much  the  optimal  shading  weight  selection  depends  on  local  noise  versus  array 
characteristics,  the  array  from  figure  8  was  extended  to  double  its  size  and  the  new  shading 
weights  were  computed.  The  results  are  shown  in  figure  9,  where  the  left  half  of  the  array 
consists  of  the  same  components  as  those  shown  in  figure  8  and  the  same  frequency  range  and 
noise  data  were  applied  in  both  cases.  While  there  are  some  elements  that  are  practically 
“zeroed-ouf  ’  in  both  cases,  the  overall  pattern  of  the  optimal  shading  weight  selection  changes 
dramatically.  This  suggests  that  selecting  the  correct  array  shading  weights  is  more  complicated 
than  just  removing  the  noisiest  elements  or  scaling  weights  inversely  with  noise  level.  Rather, 
the  complicated  patterns  of  noise  cross-correlations  impact  the  performance  of  the  array  as  a 
spatial  filter  in  a  complex  manner.  The  ability  of  the  optimization  approach  to  expose  those 
patterns  is  responsible  for  the  tremendous  performance  improvements  that  are  garnered. 


Figure  7.  Optimal  Shading  Weights  for  80-Element  Array  with  White  Noise 
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6.  CONCLUSIONS 


A  numerical  optimization  algorithm  for  determining  the  optimal  conventional  beamformer 
array  shading  weights  has  been  developed.  This  algorithm  is  based  on  maximizing  the 
beamformer  deflection  coefficient  using  knowledge  of  the  received  noise  characteristics  of  the 
individual  array  elements.  Two  algorithms  were  developed.  The  first  used  a  robust  off-the-shelf 
numerical  optimization  method  and  proved  reliable  for  small  arrays  (<100  elements),  while  the 
second  algorithm  was  a  customized  numerical  approach  that  is  specific  to  the  array  shading 
weight  optimization  problem.  The  algorithms  were  tested  against  test  array  data  and  showed 
tremendous  improvement  (usually  greater  than  20  dB)  in  signal-to-noise  ratio  relative  to  a  unity 
shaded  array. 

The  new  array  shading  weight  optimization  algorithm  may  be  used  in  a  variety  of  ways.  For 
array  designers,  the  limits  of  conventional  delay-and-sum  beamforming  can  now  be  adequately 
explored  and  compared  to  frequency-domain  adaptive  beamforming  techniques  given  the  same 
assumptions  (i.e.,  both  approaches  take  advantage  of  the  same  observed  noise  data). 

Furthermore,  the  use  of  noise  models  in  design  arrays  with  the  optimal  shading  can  point  out  the 
proper  sizing  of  an  array  to  take  advantage  of  the  spatial  complexity  of  the  noise  field,  thus 
taking  full  advantage  of  the  natural  spatial  filtering  ability  of  arrays.  Also,  the  algorithms  could 
provide  guidance  on  the  “most  important”  sections  of  a  given  array  design  by  showing  which 
elements  receive  the  largest  optimal  weights.  Finally,  a  real-time  version  of  the  algorithms  could 
be  coded  to  provide  an  in  situ  re-shading  scheme  for  tactical  arrays  based  on  the  currently 
observed  noise  characteristics.  This  would  greatly  improve  array  performance  without  impacting 
the  overall  structure  of  the  delay-and-sum  beamformer. 
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APPENDIX 

FLOWCHART  OF  THE  ARRAY  SHADING  SQP  ALGORITHM 


This  appendix  presents  a  flowchart  of  the  AS-SQP  algorithm.  The  variables  in  the 
flowchart  are  defined  as  follows: 


N: 

M: 

W(N) : 
X(N) : 

OBJ: 

G(N) : 
H(N,N) : 
A(M,N) : 
B(M) : 
BL(N) : 
BU(N) : 
IPRINT : 
MAXITR : 
IEQ : 

MAXMAJ  : 
STEPU  : 
STEPI : 

WDIFF  : 

IMAJ : 


Number  of  unknowns  (hydrophones) 

Number  of  constraints  (one  equality  constraint) 

Hydrophone  weight  vector 

Search  direction  vector 

Objective  function 

Objective  function  gradient  vector 

Hessian  matrix  of  the  objective  function 

Constraint  matrix  (each  element  =1) 

Vector  containing  RHS  of  constraints  ( =  1) 

Vector  containing  lower  bound  of  search  direction  vector 
Vector  containing  upper  bound  of  search  direction  vector 
Output  indicator  for  subroutine  DIQP 

Maximum  number  of  iterations  permitted  in  subroutine  DIQP  ( =  3*N) 

Number  of  equality  constraints  ( =  1) 

Maximum  number  of  major  iterations  permitted  for  optimization 
User-defined  step  size 

Variable  containing  the  maximum  step  size  allowed  by  weight  vector  at  current 
iterate 

Difference  between  current  and  previous  iterates  of  weight  vector 
(sum  of  absolute  values  of  components  of  difference  vector) 

Major  iteration  counter 


The  program  starts  by  reading  the  number  of  unknowns  and  constraints.  The  lower-  and  upper- 
bound  constraints  for  the  components  of  the  search  direction  vector  are  initially  set  to  -1  and  1, 
respectively,  in  accordance  with  condition  (20c).  After  the  constraint  variables  and  program 
parameters  are  set,  the  weight  vector  and  search  direction  vector  are  initialized.  The 
optimization  commences  with  the  evaluation  of  the  objective  function,  gradient,  and  Hessian  at 
the  initial  point  w(0).  These  quantities,  derived  in  section  2,  are  implemented  in  subroutine 
FNOBJ.  After  the  objective  function,  gradient,  and  Hessian  are  evaluated,  the  QP  subproblem 
(20)  is  carried  out  through  use  of  subroutine  DIQP,  developed  by  Bunch  and  Kaufman  (reference 

7). 


After  evaluation  of  the  QP  subproblem,  the  step  size  is  determined  in  accordance  with 
formula  (33).  Upon  determination  of  the  step  size,  the  weight  vector  is  updated  via  formula  (26). 
For  each  component  of  the  weight  vector  that  is  active  (i.e.,  equal  to  zero),  the  lower  bound  for 
the  corresponding  component  of  the  search  direction  vector  is  set  to  zero.  For  the  remaining 
components,  the  lower  bound  is  set  to  -1 .  The  major  iteration  counter  IMAJ  is  then  updated,  and 
the  sum  of  the  absolute  values  of  the  differences  of  the  weight  vector  components  at  the  current 
and  previous  iterates  is  evaluated  and  stored  in  the  variable  WDIFF;  i.e., 


A-l 


M 

WDIFF  =  £  | w'f+1)  -  w™  | .  (A-1) 

m= 1 

If  WDIFF  is  less  than  a  preset  tolerance  for  two  successive  iterations,  an  optimum  solution  is 
attained  and  the  program  terminates.  In  addition,  if  the  maximum  number  of  iterations  is 
achieved  before  an  optimum  solution  is  obtained  (i.e.,  IMAJ  =  MAXMAJ),  the  program 
terminates.  The  above  procedure  is  repeated  until  either  a  convergent  result  is  obtained  or  the 
maximum  allowed  number  of  iterations  is  reached.  The  computer  program  is  written  in 
FORTRAN  95  and  is  run  in  double  precision. 
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